\documentclass{article}

\usepackage{spconf}

% Title.

\title{Multimodal Brain Tumor Segmentation}

% Authors and Single addres

\name{Marco Zanger, Jorge Calvar and Demian Wasserman}

%\thanks{This work was supported by...}

\address{Common address, department \\
         City, etc \\
         optional e-mail address}

% ---------- Definicion de elementos
\def\E{\mathcal{E}}
\def\B{\mathcal{B}}
\def\W{\mathcal{W}}
\def\K{\mathcal{K}}
\def\R{\mathcal{R}}
\def\P{\mathcal{P}}
\def\D{\mathcal{D}}
\def\N{\mathcal{N}}
\def\G{\mathcal{G}}
\def\V{\mathcal{V}}
\def\L{\mathcal{L}}
\def\S{\mathcal{S}}
\def\R{\mathcal{R}}
\def\O{\mathcal{O}}

\begin{document}
	
%\ninept
%
\maketitle
%

\begin{abstract}
Automating MRI data segmentation is challenging due to the unknown characterization for tumor and edema. In this paper we propose an approach based in energy minimization focused in using different modalities (T1, T2, Flair, etc) at the same time. Extending previous works to make use of all the information stored not only in one image we are able to extend precision for region characterization. Tests are realized comparing results against manually segmented tumors made by experts.
\end{abstract}

\begin{keywords}
Segmentation, Tumour, Multimodal, \\ Energy Minimization, MRI
\end{keywords}

\section{Introduction}
\label{sec:intro}

Brain tumour and edema segmentation from MRI data is a time consuming and complex task only performed manually by experts of the field and very important in patients diagnosis. Automating this process and obtaining measurable results is hard due to the high diversity of tumors and their properties; tissue size, disposition or shape are just some of many aspects that vary from case to case resulting in an unknown characterization and, some times, hard to difference from normal tissues. Some benefits from introducing automatic processing in this field are: to reduce dependency between results and labelers, reproduce analysis very easily, being able to keep track of patients status through time and to work with 3D images and multiple modalities without increasing time results and complexity. \\
One of the most used techniques in image segmentation is the energy minimization. This metric comes from the need to compare one segmentation to another, allowing to measure how well the segmentation approaches the input image. Other  models applied in this area are: snake \cite{Kass1988}, active contours \cite{Isard1998} and shortest path technics \cite{Mortensen1998}.\\
The model we use in this work is based on \cite{Mumford1989}. Using a simple probabilistic function implemented over graph cuts and mac-flow technics, we make use of the information stored in different modalities in order to be more accurate characterizing regions. This results in an easy to understand and implement algorithm based in a graph cut technic which is known to be correct and efficient. Extending the number of input images to characterize regions of interest we avoid the problem that comes from trying to identify different regions that do not differ in intensities in only one image (in T1 with contrast liquid and blood cannot be differenced). Studies like \cite{Cobzas2007} present an algorithm based not only on intensities due to the lack of information to distinguish some tissues from others, but also on texture. In our approach we assume that MRI images do behave, in most cases, near to a peace-wise smooth image with no texture for almost all tissues.

\section{Problem formulation}

Our task consist in choosing the best from all possible image approximations $f$ using the well known energy expression

\begin{equation} \label{eq:energy}
E(f) = \sum_{\{p,q\} \in N} V_{p,q}(f_p,f_q) + \sum_{p \in P}D_p(f_p)
\end{equation}

where $V$ is the penalty function working on our $N$ set of interacting voxels in charge of measuring $f$ approximation to a piecewise smooth function, while $D_p$ measures how well the label $f_p$ fits the voxel $p$ from the input image. Then our task consist in finding a regular region partition $\Omega = \{\Omega_i | i=1 \dots n\}$ that represents our regular regions of interest. Our first choice was to assume each voxels intensity as a result of a well known random process. Taking $p_i$ as the probability for a voxel to belong to region $\Omega_i$ and following \cite{Paragios2002}, we know that maximizing each voxels probability is the same as minimizing the energy

\begin{equation} \label{eq:probabilidad}
E(\Omega, \theta) = - \sum_i
\int_{\Omega_i} \log p_i(I(x)|\theta_i) dx + \nu|C|
\end{equation}

In order to obtain expression (\ref{eq:probabilidad}) some simplifications had been made, and these lead
to certain restrictions over the model. First, the number of regions to identify must be 
known from the very beginning, preventing us from creating an unattended procedure. Also as voxel intensities within the same region are supposed to be the result of the same random process, identically distributed and independent from each other, our model expects input images to be as close as possible to a piecewise smooth image.

As we said, our main simplification is to assume that every region intensity
distribution comes from the same density function but with different parameters $\theta = \{\theta_i, \ldots,\theta_N\}$. Knowing the amounts of regions to identify and the parameters for each density
function we can express our energy minimization as the following

\begin{equation}
E(\Omega, \theta) = - \sum_i
\int_{\Omega_i} \log p_i(I(x)|\theta_i) dx + \nu|C|
\end{equation}

where $p(I|\theta_i)$ is the a posteriori probability for region $\Omega_i$ with
parameter $\theta_i$. One of the simplest approaches in order to characterize
these parameters and work with different modalities at the same time is to work with multivariate Gaussian density as an approximation of region statistics \cite{Rousson2004}:

\begin{equation}
p(z|\mu, \Sigma) =
\frac{1}{(2\pi)^{d/2}|\Sigma|}e^{-\frac{1}{2}(z-\mu)^T\Sigma^{-1}(z-\mu)}
\end{equation}

in \cite{fischl2002} we find a voxels intensities distribution analysis for almost all tissues obtained from a MRI study, showing most of them behave near a normal distribution or as a sum of several normal distribution.

\subsection{Energy Minimization}

In order to minimize our energy function iteratively, we decided to use an Expectation Minimization algorithm. This consist in defining two steps; the first one, step E, fixes the probabilistic parameters to find the best possible labeling, then the second step M uses the labeling obtained in E to find the best probabilistic parameters.

\subsubsection{Finding the best labeling}

This step is addressed using one of the well known fast approximate energy minimization algorithms based in graph cuts. \cite{Boykov2001} describes two different methods for obtaining local minimum and the properties that results from using one or the other. These different methods are based in what they define as moves, which represents the change of one or more voxels from one region to another. Within these possible methods there is one called expansion move, which is described to be the most effective one assuming one restriction over the model. The technique they describe solves the problem assuming that (each) $V$ is a metric and, thus, satisfies the triangle inequality \cite{Boykov2001}.
The first step of our method consist in adding initial points corresponding to the different regions to identify. In order to do this the user selects samples of the image that characterize each region (view image ??) %\ref{im:initipoints})

%\begin{figure}[h]
%  \centering
%  \includegraphics[width=100px]{}
%  \caption{Initial Points} \label{im:initpoints}
%\end{figure}

using these samples we are able to obtain the function density parameters to assign probabilities for each voxel intensity to belong to a region. The next step consist in building the graph following \cite{Boykov2001}. First, we define our neighborhood configuration to include all the direct connected voxels. This result in an up to twenty six neighbors for each voxel. 

\subsubsection{Obtaining new parameters}

After fixing the labeling $f$ we obtain new parameters $\theta$. As we defined using the normal distribution for all the regions, this step is fairly simple, and consists of computing the parameters for each region.

\section{Experimental results}

The output of our experiments consist in a set of spatially coherent regions. In order to measure the accurate of our method, we ask for an expert to segment the tumor by hand for a set of images and compare these segmentation against our results. To do this we used the Dice coefficient, which is a similarity measure, in function of the smooth factor and the number of iterations.

\subsection{Initialization}
After several tests we found out that the initialization has a significant effect over the final results. 

\section{Conclusions and Future Work}

We have presented a new technique for segmentation based in the utilization of different modalities and the information within. As we expected, our main assumption over tumors intensity distribution is the weakest of all. As these tissues do not present a known characterization one of the optimization over our model would be to use a Student t-distribution. This distribution arises when (as in nearly all practical statistical work) the population standard deviation is unknown and has to be estimated from the data, or there are a certain number of outsiders that must be taken into account.



\bibliographystyle{IEEEbib}
\bibliography{paper}

\end{document}